Set Up¶
# Import Libraries
import json # Reading and writing JSON
from glob import glob # File pattern matching
import os # Operating system access
import pathlib # File navigation
import holoviews as hv # Interactable plot
import panel as pn # Plot formatting
pn.extension() # Plot formatting
hv.extension('bokeh') # Interactable plot
import numpy as np # Numerical computing
import earthpy # Work with geospatial data
import earthpy.api.appeears as eaapp # Access earthpy api
import xarray as xr # Multi-dimensional arrays
import hvplot.xarray # HVPlot and xarray
import rioxarray as rxr # Raster for xarray
import matplotlib.pyplot as plt # Plotting
import statsmodels.formula.api as smf # Statistical models
import pandas as pd # Work with tabular data
import geopandas as gpd # Wowk with geospatial data
import hvplot.pandas # Plot geospatial data
# Establish project directory
project = earthpy.Project(
'Gila River Vegetation', dirname='NDVI_data')
project.get_data()
Establish Boundaries¶
# Load in the boundary data
aitsn_gdf = gpd.read_file(project.project_dir / 'tl_2020_us_aitsn')
print(type(aitsn_gdf))
aitsn_gdf.head()
<class 'geopandas.geodataframe.GeoDataFrame'>
| AIANNHCE | TRSUBCE | TRSUBNS | GEOID | NAME | NAMELSAD | LSAD | CLASSFP | MTFCC | FUNCSTAT | ALAND | AWATER | INTPTLAT | INTPTLON | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 2430 | 653 | 02419073 | 2430653 | Red Valley | Red Valley Chapter | T2 | D7 | G2300 | A | 922036695 | 195247 | +36.6294607 | -109.0550394 | POLYGON ((-109.2827 36.64644, -109.28181 36.65... |
| 1 | 2430 | 665 | 02419077 | 2430665 | Rock Point | Rock Point Chapter | T2 | D7 | G2300 | A | 720360268 | 88806 | +36.6598701 | -109.6166836 | POLYGON ((-109.85922 36.49859, -109.85521 36.5... |
| 2 | 2430 | 675 | 02419081 | 2430675 | Rough Rock | Rough Rock Chapter | T2 | D7 | G2300 | A | 364475668 | 216144 | +36.3976971 | -109.7695183 | POLYGON ((-109.93053 36.40672, -109.92923 36.4... |
| 3 | 2430 | 325 | 02418975 | 2430325 | Indian Wells | Indian Wells Chapter | T2 | D7 | G2300 | A | 717835323 | 133795 | +35.3248534 | -110.0855000 | POLYGON ((-110.24222 35.36327, -110.24215 35.3... |
| 4 | 2430 | 355 | 02418983 | 2430355 | Kayenta | Kayenta Chapter | T2 | D7 | G2300 | A | 1419241065 | 1982848 | +36.6884391 | -110.3045616 | POLYGON ((-110.56817 36.73489, -110.56603 36.7... |
# Plot GRIC and SRIC boundary overview
overview_plot = aitsn_gdf.hvplot(
geo=True, tiles='CartoLight',
frame_width=500,
legend=False, fill_color=None, fill_alpha = .5,
line_color='blue',
colorbar = False,
hover_cols='all',
xlim=(-112.4, -111.4),
ylim=(32.9, 33.6),
title = 'GRIC & SRIC Boundary Overview'
)
overview_plot
# Save the plot as an HTML file
hv.save(overview_plot, 'img/overview_plot.html', fmt='html')
WARNING:bokeh.core.validation.check:W-1005 (FIXED_SIZING_MODE): 'fixed' sizing mode requires width and height to be set: figure(id='44422e10-b12c-4191-802e-c44f25f55e8f', ...)
# Plot GRIC Boundary and set up gdf
gric_gdf = aitsn_gdf.loc[aitsn_gdf.AIANNHCE=='1310'].dissolve()
gric_bound_plot = gric_gdf.hvplot(
geo=True, tiles='CartoLight',
fill_color='green', line_color='Blue',
fill_alpha=.3, # must be b/w 0-1
title='Gila River Indian Community',
frame_width=500)
gric_bound_plot
# Plot SRIC Boundary and set up gdf
sric_gdf = aitsn_gdf.loc[aitsn_gdf.AIANNHCE=='3340'].dissolve()
sric_bound_plot = sric_gdf.hvplot(
geo=True, tiles='CartoLight',
fill_color='green', line_color='Blue',
fill_alpha=.3, # must be b/w 0-1
title='Salt River Indian Community',
frame_width=500)
sric_bound_plot
# Save both boundary plots
hv.save(gric_bound_plot, 'img/gric_bound_plot.html', fmt='html')
hv.save(sric_bound_plot, 'img/sric_bound_plot.html', fmt='html')
WARNING:bokeh.core.validation.check:W-1005 (FIXED_SIZING_MODE): 'fixed' sizing mode requires width and height to be set: figure(id='ca7dbb49-e1f2-452c-b4ba-b2ab2bd90c48', ...) WARNING:bokeh.core.validation.check:W-1005 (FIXED_SIZING_MODE): 'fixed' sizing mode requires width and height to be set: figure(id='0f6af472-c0d0-4107-abca-dd4d2e9b0375', ...)
Download NDVI Data¶
The GRIC NDVI data is included in the download for the AITSN Boundary data used above, and the earthpy API can be used to download the SRIC NDVI data using the NASA Earthdata API
# Set up download key
ndvi_downloader = eaapp.AppeearsDownloader(
### give your download a name
download_key = "salt-river-ndvi",
### tell it to put the data in your project that you already defined
project = project,
### specify the MODIS product you want
product = 'MOD13Q1.061',
layer = '_250m_16_days_NDVI',
### choose a start date and end data
start_date = "05-24",
end_date = "08-29",
### recurring means you want those dates over multiple years
recurring = True,
### specify the range of years you want
year_range = [2001, 2022],
### specify the polygon you want to get NDVI data for
polygon = sric_gdf
)
# Download files if not already saved
ndvi_downloader.download_files(cache=True)
<generator object Path.rglob at 0x788728006bd0>
# Get a sorted list of file path for only the 'gila_river-ndvi' folder
gila_dir = project.project_dir / 'gila-river-ndvi'
gric_paths = sorted(list(gila_dir.rglob('*NDVI*.tif')))
# Display the first and last three files paths to check the pattern
gric_paths[:3], gric_paths[-3:]
([PosixPath('/workspaces/data/ndvi_data/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2001145000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_data/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2001161000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_data/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2001177000000_aid0001.tif')],
[PosixPath('/workspaces/data/ndvi_data/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2022209000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_data/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2022225000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_data/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2022241000000_aid0001.tif')])
# Get a sorted list of file path for only the 'satl-river-ndvi' folder
salt_dir = project.project_dir / 'salt-river-ndvi'
sric_paths = sorted(list(salt_dir.rglob('*NDVI*.tif')))
# Display the first and last three files paths to check the pattern
sric_paths[:3], sric_paths[-3:]
([PosixPath('/workspaces/data/ndvi_data/salt-river-ndvi/MOD13Q1.061_2001129_to_2022241/MOD13Q1.061__250m_16_days_NDVI_doy2001129000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_data/salt-river-ndvi/MOD13Q1.061_2001129_to_2022241/MOD13Q1.061__250m_16_days_NDVI_doy2001145000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_data/salt-river-ndvi/MOD13Q1.061_2001129_to_2022241/MOD13Q1.061__250m_16_days_NDVI_doy2001161000000_aid0001.tif')],
[PosixPath('/workspaces/data/ndvi_data/salt-river-ndvi/MOD13Q1.061_2001129_to_2022241/MOD13Q1.061__250m_16_days_NDVI_doy2022209000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_data/salt-river-ndvi/MOD13Q1.061_2001129_to_2022241/MOD13Q1.061__250m_16_days_NDVI_doy2022225000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_data/salt-river-ndvi/MOD13Q1.061_2001129_to_2022241/MOD13Q1.061__250m_16_days_NDVI_doy2022241000000_aid0001.tif')])
# Compare datasets
print(len(sric_paths))
print(len(gric_paths))
171 154
The Salt River NDVI files included an extra day on the front end of some years (day 119). We will filter both file lists by the same date range later to ensure everything lines up.
Format NDVI Data¶
# Add a date dimension to each file from the file name
doy_start = -25
doy_end = -19
all_paths = {
"gric": gric_paths,
"sric": sric_paths
}
results = {
"gric": [],
"sric": []
}
for label, paths in all_paths.items():
for path in paths:
# Get date from file name
doy = path.name[doy_start:doy_end]
date = pd.to_datetime(doy, format='%Y%j')
# Open dataset
da = rxr.open_rasterio(path, mask_and_scale=True).squeeze()
# Add date dimension and clean metadata
da = da.assign_coords({'date': date})
da = da.expand_dims({'date': 1})
da.name = 'NDVI'
# Append to the correct list
results[label].append(da)
gric_das = results["gric"]
sric_das = results["sric"]
print(len(gric_das))
print(len(sric_das))
154 171
# Combine NDVI images from all dates
combined = {}
for label, da_list in results.items():
combined[label] = xr.combine_by_coords(
da_list, coords=['date'], compat='override'
)
# Access the final outputs
gric_da = combined["gric"]
sric_da = combined["sric"]
print(sric_da)
print(gric_da)
<xarray.Dataset> Size: 6MB
Dimensions: (date: 171, y: 71, x: 130)
Coordinates:
band int64 8B 1
* x (x) float64 1kB -111.9 -111.9 -111.9 ... -111.6 -111.6 -111.6
* y (y) float64 568B 33.58 33.58 33.58 33.58 ... 33.44 33.44 33.44
spatial_ref int64 8B 0
* date (date) datetime64[ns] 1kB 2001-01-12 2001-01-14 ... 2022-01-24
Data variables:
NDVI (date, y, x) float32 6MB 0.2532 0.2532 0.2532 ... 0.18 0.1789
<xarray.Dataset> Size: 48MB
Dimensions: (date: 154, y: 203, x: 382)
Coordinates:
band int64 8B 1
* x (x) float64 3kB -112.3 -112.3 -112.3 ... -111.5 -111.5 -111.5
* y (y) float64 2kB 33.39 33.39 33.38 33.38 ... 32.97 32.97 32.97
spatial_ref int64 8B 0
* date (date) datetime64[ns] 1kB 2001-01-14 2001-01-16 ... 2022-01-24
Data variables:
NDVI (date, y, x) float32 48MB 0.8282 0.6146 ... 0.2146 0.2085
# Make sure each data array only contains data for dates present in both
common_dates = np.intersect1d(gric_da.date.values, sric_da.date.values)
gric_matched = gric_da.sel(date=common_dates)
sric_matched = sric_da.sel(date=common_dates)
len(gric_matched.date) - len(sric_matched.date)
0
# Clip the data arrays to only include data within each boundary
gric_ndvi = gric_matched.rio.clip(gric_gdf.geometry, from_disk=True)
sric_ndvi = sric_matched.rio.clip(sric_gdf.geometry, from_disk=True)
# Calculate spatial means for each dataset
gric_mean = gric_ndvi.mean(dim=['x','y'])
sric_mean = sric_ndvi.mean(dim=['x','y'])
# Check everything looks as expected
gric_mean
<xarray.Dataset> Size: 2kB
Dimensions: (date: 154)
Coordinates:
band int64 8B 1
* date (date) datetime64[ns] 1kB 2001-01-14 2001-01-16 ... 2022-01-24
spatial_ref int64 8B 0
Data variables:
NDVI (date) float32 616B 0.1927 0.1838 0.1934 ... 0.2401 0.2125Set Up DiD Analysis¶
A difference in difference OLS regression assumes we have similar trends across each dataset, which is largely true here due to proximity but is roughly tested below. Once similar trends can be established, effects from outside of shared trends should be easier to identify. More information on DiD can be found here.
# Modify arrays to dataframes and combine
gric_df = gric_mean['NDVI'].to_dataframe().reset_index()
sric_df = sric_mean['NDVI'].to_dataframe().reset_index()
gric_df['group'] = 'gric'
sric_df['group'] = 'sric'
df = pd.concat([gric_df, sric_df], ignore_index=True)
df['year'] = df['date'].dt.year
df['month'] = df['date'].dt.month
df['doy'] = df['date'].dt.dayofyear
df['treated'] = (df['group'] == 'gric').astype(int)
df['post'] = (df['year'] >= 2005).astype(int)
df['did'] = df['treated'] * df['post']
df
| date | band | spatial_ref | NDVI | group | year | month | doy | treated | post | did | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 2001-01-14 | 1 | 0 | 0.192719 | gric | 2001 | 1 | 14 | 1 | 0 | 0 |
| 1 | 2001-01-16 | 1 | 0 | 0.183815 | gric | 2001 | 1 | 16 | 1 | 0 | 0 |
| 2 | 2001-01-17 | 1 | 0 | 0.193379 | gric | 2001 | 1 | 17 | 1 | 0 | 0 |
| 3 | 2001-01-19 | 1 | 0 | 0.206412 | gric | 2001 | 1 | 19 | 1 | 0 | 0 |
| 4 | 2001-01-20 | 1 | 0 | 0.205276 | gric | 2001 | 1 | 20 | 1 | 0 | 0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 303 | 2022-01-17 | 1 | 0 | 0.203360 | sric | 2022 | 1 | 17 | 0 | 1 | 0 |
| 304 | 2022-01-19 | 1 | 0 | 0.194895 | sric | 2022 | 1 | 19 | 0 | 1 | 0 |
| 305 | 2022-01-20 | 1 | 0 | 0.208739 | sric | 2022 | 1 | 20 | 0 | 1 | 0 |
| 306 | 2022-01-22 | 1 | 0 | 0.245338 | sric | 2022 | 1 | 22 | 0 | 1 | 0 |
| 307 | 2022-01-24 | 1 | 0 | 0.238629 | sric | 2022 | 1 | 24 | 0 | 1 | 0 |
308 rows × 11 columns
# Clean up dataframe
df = df.drop(columns=['band', 'spatial_ref'])
df = df.rename(columns={'NDVI': 'ndvi'})
df = df.reset_index(drop=True)
df
| date | ndvi | group | year | month | doy | treated | post | did | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 2001-01-14 | 0.192719 | gric | 2001 | 1 | 14 | 1 | 0 | 0 |
| 1 | 2001-01-16 | 0.183815 | gric | 2001 | 1 | 16 | 1 | 0 | 0 |
| 2 | 2001-01-17 | 0.193379 | gric | 2001 | 1 | 17 | 1 | 0 | 0 |
| 3 | 2001-01-19 | 0.206412 | gric | 2001 | 1 | 19 | 1 | 0 | 0 |
| 4 | 2001-01-20 | 0.205276 | gric | 2001 | 1 | 20 | 1 | 0 | 0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 303 | 2022-01-17 | 0.203360 | sric | 2022 | 1 | 17 | 0 | 1 | 0 |
| 304 | 2022-01-19 | 0.194895 | sric | 2022 | 1 | 19 | 0 | 1 | 0 |
| 305 | 2022-01-20 | 0.208739 | sric | 2022 | 1 | 20 | 0 | 1 | 0 |
| 306 | 2022-01-22 | 0.245338 | sric | 2022 | 1 | 22 | 0 | 1 | 0 |
| 307 | 2022-01-24 | 0.238629 | sric | 2022 | 1 | 24 | 0 | 1 | 0 |
308 rows × 9 columns
# Plot mean yearly NDVI across both to compare trends prior to 2005
df_yearly = df.groupby(['year', 'group']).ndvi.mean().reset_index()
plt.figure(figsize=(10,6))
for g in ['gric', 'sric']:
subset = df_yearly[df_yearly['group'] == g]
plt.plot(subset['year'], subset['ndvi'], marker='o', label=g.upper())
# vertical line at 2005
plt.axvline(2005, color='k', linestyle='--', label='Post Settlement')
plt.title("Parallel Trends Check: NDVI on GRIC compared to SRIC")
plt.ylabel("Mean NDVI")
plt.xlabel("Year")
plt.legend()
plt.grid(True)
plt.savefig("img/ndvi_parallel_trends.png", dpi=300, bbox_inches="tight")
plt.show()
# Run DID model
did_model = smf.ols(
formula="ndvi ~ treated*post + C(year) + C(month)",
data=df
).fit()
print(did_model.summary())
OLS Regression Results
==============================================================================
Dep. Variable: ndvi R-squared: 0.584
Model: OLS Adj. R-squared: 0.550
Method: Least Squares F-statistic: 17.33
Date: Wed, 10 Dec 2025 Prob (F-statistic): 1.55e-41
Time: 20:36:14 Log-Likelihood: 805.33
No. Observations: 308 AIC: -1563.
Df Residuals: 284 BIC: -1473.
Df Model: 23
Covariance Type: nonrobust
===================================================================================
coef std err t P>|t| [0.025 0.975]
-----------------------------------------------------------------------------------
Intercept 0.2310 0.006 41.920 0.000 0.220 0.242
C(year)[T.2002] -0.0213 0.007 -3.052 0.002 -0.035 -0.008
C(year)[T.2003] -0.0055 0.007 -0.785 0.433 -0.019 0.008
C(year)[T.2004] -0.0199 0.007 -2.855 0.005 -0.034 -0.006
C(year)[T.2005] 0.0385 0.005 8.022 0.000 0.029 0.048
C(year)[T.2006] 0.0120 0.005 2.491 0.013 0.003 0.021
C(year)[T.2007] -0.0158 0.005 -3.298 0.001 -0.025 -0.006
C(year)[T.2008] 0.0141 0.005 2.941 0.004 0.005 0.024
C(year)[T.2009] -0.0116 0.005 -2.417 0.016 -0.021 -0.002
C(year)[T.2010] 0.0104 0.005 2.162 0.031 0.001 0.020
C(year)[T.2011] -0.0207 0.005 -4.319 0.000 -0.030 -0.011
C(year)[T.2012] -0.0211 0.005 -4.390 0.000 -0.031 -0.012
C(year)[T.2013] -0.0063 0.005 -1.310 0.191 -0.016 0.003
C(year)[T.2014] -0.0022 0.005 -0.458 0.647 -0.012 0.007
C(year)[T.2015] -0.0025 0.005 -0.517 0.606 -0.012 0.007
C(year)[T.2016] -0.0183 0.005 -3.822 0.000 -0.028 -0.009
C(year)[T.2017] -0.0036 0.005 -0.755 0.451 -0.013 0.006
C(year)[T.2018] -0.0241 0.005 -5.030 0.000 -0.034 -0.015
C(year)[T.2019] 0.0085 0.005 1.780 0.076 -0.001 0.018
C(year)[T.2020] 0.0101 0.005 2.094 0.037 0.001 0.019
C(year)[T.2021] 0.0271 0.005 5.650 0.000 0.018 0.037
C(year)[T.2022] -0.0029 0.005 -0.608 0.544 -0.012 0.007
treated -0.0341 0.005 -6.916 0.000 -0.044 -0.024
post -0.0086 0.005 -1.572 0.117 -0.019 0.002
treated:post 0.0085 0.005 1.561 0.120 -0.002 0.019
==============================================================================
Omnibus: 190.755 Durbin-Watson: 1.114
Prob(Omnibus): 0.000 Jarque-Bera (JB): 3348.597
Skew: 2.173 Prob(JB): 0.00
Kurtosis: 18.558 Cond. No. 5.76e+15
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The smallest eigenvalue is 2.14e-29. This might indicate that there are
strong multicollinearity problems or that the design matrix is singular.
# Save the model output
with open("img/did_model_summary.html", "w") as f:
f.write(did_model.summary().as_html())
DiD Conclusion¶
Unfortunately, this type of analysis does not lead to very strong results, largely due to annual variability overshadowing any effects the increased flow might have had. Our primary variable is treated:post = 0.0085, which on a scale of 0-1 is a very small (although still positive!) impact when taking out similar trends between the datasets. I believe this is also partly due to the GRIC seeing increased development as time goes on. Having more data on the front end could help as well. By analyzing negative growth trends compared to street maps or similar information, my next step would be to mask over any increased development in both the GRIC and SRIC to isolate as much 'pure' vegetation data for analysis as possible, but I have not been able to work that out yet beyond estimation and eyeballing maps.¶
# Plot the two areas comparing mean values before and after 2005 for overview
# GRIC NDVI
gric_pre = gric_ndvi.sel(date=slice("2001-01-01", "2004-12-31")).NDVI.mean(dim='date')
gric_post = gric_ndvi.sel(date=slice("2005-01-01", "2022-12-31")).NDVI.mean(dim='date')
gric_diff = gric_post - gric_pre
# SRIC NDVI (control)
sric_pre = sric_ndvi.sel(date=slice("2001-01-01", "2004-12-31")).NDVI.mean(dim='date')
sric_post = sric_ndvi.sel(date=slice("2005-01-01", "2022-12-31")).NDVI.mean(dim='date')
sric_diff = sric_post - sric_pre
# Plot GRIC pre/post NDVI change
plt.figure(figsize=(5,4))
gric_diff.plot(cmap='RdYlGn', vmin=-0.05, vmax=0.05)
plt.title('GRIC NDVI change post-2005 (mean difference)')
plt.savefig("img/gric_change.png", dpi=300, bbox_inches="tight")
plt.show()
# Plot SRIC pre/post NDVI change
plt.figure(figsize=(5,4))
sric_diff.plot(cmap='RdYlGn', vmin=-0.05, vmax=0.05)
plt.title('SRIC NDVI change post-2005 (mean difference)')
plt.savefig("img/sric_change.png", dpi=300, bbox_inches="tight")
plt.show()
Plot Running Average Against Pre-Settlement Mean¶
I have not figured out how to implement SRIC data as a control on a raster level, but plottting each years running average from 2005 compared to mean values prior can be helpful in identifying areas of increased and decreased NDVI data that would be worth looking into further.
# Calculate mean values prior to the settlement
gric_ndvi_pre = (gric_ndvi
.sel(date=slice('2001', '2004'))
.mean('date')
.NDVI
)
# Calculate annual mean values post settlement
gric_ndvi_post = gric_ndvi.sel(date=slice('2005', '2022')).NDVI
gric_post_yearly = gric_ndvi_post.groupby('date.year').mean('date')
# Confirm grouping correctly
len(gric_post_yearly)
18
# Calculate difference between each year post and mean prior
gric_diff = gric_post_yearly - gric_ndvi_pre
print(gric_diff.dims)
print(gric_diff.coords)
('year', 'y', 'x')
Coordinates:
band int64 8B 1
* x (x) float64 3kB -112.3 -112.3 -112.3 ... -111.5 -111.5 -111.5
* y (y) float64 2kB 33.39 33.38 33.38 33.38 ... 32.97 32.97 32.97
spatial_ref int64 8B 0
* year (year) int64 144B 2005 2006 2007 2008 ... 2019 2020 2021 2022
# Plot the difference
diff_plot = (
gric_diff.hvplot(x='x', y='y', cmap='RdYlGn', geo=True, groupby='year',
clim=(-0.5, 0.5),
title = 'Gila River Indian Community NDVI\n'
'Comparing 2001-2004 to Following Years')
*
gric_gdf.hvplot(geo=True, fill_color=None, line_color='black')
)
diff_plot
BokehModel(combine_events=True, render_bundle={'docs_json': {'e9c377f4-bb84-4fea-81f6-831fb75b2eac': {'version…
# Move the slider to the bottom of the plot and force 2005 start
years = list(diff_plot.kdims[0].values)
years.sort() # just to be safe
diff_plot_bottom = pn.panel(
diff_plot,
widget_location='bottom',
widgets={
'year': pn.widgets.DiscreteSlider(
name='year',
options=years,
value=years[0]
)
}
)
diff_plot_bottom
BokehModel(combine_events=True, render_bundle={'docs_json': {'ecbbcbf4-6c17-4aca-b2db-31f9c4152d52': {'version…
# Save plot for comparison
diff_plot_bottom.save('img/gric_diff_plot.html', embed=True)
WARNING:bokeh.core.validation.check:W-1005 (FIXED_SIZING_MODE): 'fixed' sizing mode requires width and height to be set: figure(id='4137e24d-9ed0-4aa1-ac44-5883760b4de1', ...)
# Calculate cumulative sum for each year from settlement
cumulative_sum = gric_post_yearly.cumsum('year')
len(cumulative_sum)
18
# Calculate year count for each year post settlement
year_count = xr.DataArray(
range(1, len(gric_post_yearly.year) + 1),
coords=[gric_post_yearly.year],
dims=['year']
)
len(year_count)
18
# Calculate average for each year
gric_run_avg = cumulative_sum / year_count
len(gric_run_avg)
18
# Calculate average difference from pre-settlement mean
new_gric_diff = gric_run_avg - gric_ndvi_pre
print(new_gric_diff.dims)
print(new_gric_diff.coords)
('year', 'y', 'x')
Coordinates:
band int64 8B 1
* x (x) float64 3kB -112.3 -112.3 -112.3 ... -111.5 -111.5 -111.5
* y (y) float64 2kB 33.39 33.38 33.38 33.38 ... 32.97 32.97 32.97
spatial_ref int64 8B 0
* year (year) int64 144B 2005 2006 2007 2008 ... 2019 2020 2021 2022
# Plot the new running average
new_diff_plot = (
new_gric_diff.hvplot(x='x', y='y', cmap='RdYlGn', geo=True, groupby='year',
clim=(-0.5, 0.5),
title = 'Gila River Indian Community NDVI\n'
'Comparing 2001-2004 to Average of Following Years')
*
gric_gdf.hvplot(geo=True, fill_color=None, line_color='black')
)
years = list(diff_plot.kdims[0].values)
years.sort()
# Move slider to the bottom and force correct start
new_diff_plot_bottom = pn.panel(
new_diff_plot,
widget_location='bottom',
widgets={
'year': pn.widgets.DiscreteSlider(
name='year',
options=years,
value=years[0]
)
}
)
new_diff_plot_bottom
BokehModel(combine_events=True, render_bundle={'docs_json': {'15792bcd-eea6-4856-b4f2-f6dff4474570': {'version…
# Save the new plot
new_diff_plot_bottom.save('img/gric_avg_diff_plot.html', embed=True)
WARNING:bokeh.core.validation.check:W-1005 (FIXED_SIZING_MODE): 'fixed' sizing mode requires width and height to be set: figure(id='d80134e9-4371-46cf-bd91-a5df779d3367', ...)